广义摆动方程吸引域边界重整化

标准形式

广义摆动方程的标准形式为：

{\begin{cases} \begin{aligned} \frac{d}{d t} x (t) = - y (t) \\ \frac{d}{d t} y (t) = \underset{平衡项}{\underset{⏟}{- β}} + \underset{非线性恢复力}{\underset{⏟}{\sin (x)}} + \underset{速度依赖阻尼}{\underset{⏟}{(α \cos (x) - γ) y (t)}} \end{aligned} \end{cases} 1.

平衡性分析

雅可比矩阵为：

\begin{array}{r} J (x, y) = [\begin{array}{c} 0 & - 1 \\ \cos (x) & α \cos (x) - γ \end{array}] \end{array}

平衡点分析

平衡点 $(x^{*}, y^{*})$ 要满足

y^{*} = 0, β - \sin x^{*} = 0 ⟹ \sin x^{*} = β .

若 $| β | < 1$ ，则有两组解：

x_{1}^{*} = \arcsin β, x_{2}^{*} = π - \arcsin β, y^{*} = 0.

注意 $\cos x_{1}^{*} = \sqrt{1 - β^{2}} > 0$ ， $\cos x_{2}^{*} = - \sqrt{1 - β^{2}} < 0$ 。
当 $\cos x^{*} < 0$ 时 $\det J (x^{*}, 0) > 0$ ，对应鞍点，因此选

(x_{s}, y_{s}) = (π - \arcsin β, 0) .

鞍点处的雅可比矩阵

把 $(x, y) = (x_{s}, 0)$ 代入 $J (x, y)$ 得

\begin{matrix} J_{s} = J (x_{s}, 0) = (\begin{matrix} 0 & - 1 \\ \sqrt{1 - β^{2}} & α \sqrt{1 - β^{2}} - γ \end{matrix}) . \end{matrix}

记迹与行列式分别为

T = tr J_{s} = α \sqrt{1 - β^{2}} - γ, D = \det J_{s} = \sqrt{1 - β^{2}} > 0.

特征值与特征向量

特征值 $λ$ 满足方程

λ^{2} - T λ + D = 0 ⟹ λ_{1, 2} = \frac{T \pm \sqrt{T^{2} - 4 D}}{2} = \frac{(α \sqrt{1 - β^{2}} - γ) \pm \sqrt{(α \sqrt{1 - β^{2}} - γ)^{2} - 4 \sqrt{1 - β^{2}}}}{2} .

由于 $D > 0$ ，故一正一负，为典型鞍点。
对应的特征向量可取

\begin{matrix} v_{i} = (\begin{matrix} 1 \\ λ_{i} \end{matrix}), i = 1, 2, \end{matrix}

稳定流形和不稳定流形

归一化后，两个特征向量取

\begin{matrix} v_{i} = (\begin{matrix} 1 \\ λ_{i} \end{matrix}), ∥ v_{i} ∥ = \sqrt{1 + λ_{i}^{2}}, {\hat{v}}_{i} = \frac{v_{i}}{∥ v_{i} ∥} . \end{matrix}

带入特征值 $λ_{1} \approx 1.37052, λ_{2} \approx - 0.66896$ 得：

\begin{aligned} ∥ v_{1} ∥ & = \sqrt{1 + (1.37052)^{2}} \approx 1.697, & {\hat{v}}_{1} & = (\begin{array}{c} 1 / 1.697 \\ 1.37052 / 1.697 \end{array}) \approx (\begin{array}{c} 0.5896 \\ 0.8077 \end{array}), \\ ∥ v_{2} ∥ & = \sqrt{1 + (- 0.66896)^{2}} \approx 1.203, & {\hat{v}}_{2} & = (\begin{array}{c} 1 / 1.203 \\ - 0.66896 / 1.203 \end{array}) \approx (\begin{array}{c} 0.8313 \\ - 0.5568 \end{array}) . \end{aligned}

因 $λ_{1} > 0$ ， ${\hat{v}}_{1}$ 对应 不稳定流形（unstable manifold）方向
因 $λ_{2} < 0$ ， ${\hat{v}}_{2}$ 对应 稳定流形（stable manifold）方向

鞍点附近的泰勒展开

\begin{array}{r} f^{(5)} (x, y) = {\begin{cases} \frac{d}{d t} x (t) = - y (t) \\ \frac{d}{d t} y (t) = \underset{线性项}{\underset{⏟}{(α \sqrt{1 - β^{2}} - γ) y}} + \underset{一阶项}{\underset{⏟}{(\sqrt{1 - β^{2}} + α β y) (x - x_{s})}} + \underset{二阶项}{\underset{⏟}{\frac{1}{2} (β - α y \sqrt{1 - β^{2}}) (x - x_{s})^{2}}} \\ + \underset{三阶项}{\underset{⏟}{\frac{1}{6} (- \sqrt{1 - β^{2}} - α β y) (x - x_{s})^{3}}} + \underset{四阶项}{\underset{⏟}{\frac{1}{24} (- β + α y \sqrt{1 - β^{2}}) (x - x_{s})^{4}}} + \underset{五阶项}{\underset{⏟}{\frac{1}{120} (\sqrt{1 - β^{2}} + α β y) (x - x_{s})^{5}}} \end{cases} \end{array}

重整化

1.坐标变换

将 $[x, y]$ 坐标变换为 $[z, w]$ 坐标，其中 $z, w$ 分别为鞍点的不稳定方向和稳定方向：

\begin{matrix} a_{1}, a_{2} = v_{1} \\ b_{1}, b_{2} = v_{2} \\ x = (a_{1} z + b_{1} w) + x_{s} \\ y = (a_{2} z + b_{2} w) + y_{s} \end{matrix}

新的方程表示为

\begin{matrix} \frac{d}{d t} [\begin{matrix} z \\ w \end{matrix}] = [\begin{matrix} 0.214072037938063 w + 0.214072037938063 z - 0.0267962588229219 {(0.506615170878591 w - z)}^{5} \cdot (0.42 w + 0.42 z + 0.8) - 0.0852959198738089 {(0.506615170878591 w - z)}^{4} \cdot (0.56 w + 0.56 z - 0.6) + 0.217206259879704 {(0.506615170878591 w - z)}^{3} \cdot (0.42 w + 0.42 z + 0.8) - 0.414837187412316 {(0.506615170878591 w - z)}^{2} (- 0.56 w - 0.56 z + 0.6) - 0.336260500140295 \cdot (0.795781982453887 w - 1.57078198245389 z) (0.42 w + 0.42 z + 0.8) \\ - 0.834072037938063 w - 0.834072037938063 z - 0.0528927287677762 {(0.506615170878591 w - z)}^{5} \cdot (0.42 w + 0.42 z + 0.8) - 0.168364322224867 {(0.506615170878591 w - z)}^{4} \cdot (0.56 w + 0.56 z - 0.6) + 0.42874014116675 {(0.506615170878591 w - z)}^{3} \cdot (0.42 w + 0.42 z + 0.8) - 0.818840830788565 {(0.506615170878591 w - z)}^{2} (- 0.56 w - 0.56 z + 0.6) - 0.663739499859705 \cdot (0.795781982453887 w - 1.57078198245389 z) (0.42 w + 0.42 z + 0.8) \end{matrix}] \end{matrix}

2.摄动展开

执行摄动展开:

\begin{matrix} z = z_{0} + ϵ z_{1} + ϵ^{1} z_{2} + ϵ^{2} z_{3} + ϵ^{3} z_{4} + ϵ^{4} z_{5} \\ w = ϵ w_{1} + ϵ^{1} w_{2} + ϵ^{2} w_{3} + ϵ^{3} w_{4} + ϵ^{4} w_{5} \end{matrix}

带入上式方程，对应阶数得到对应阶数方程，形如

= 1 阶方程 = dz/dt 方程: Eq(Derivative(z1(t), t), 0.636625585963109z1(t)) dw/dt 方程: Eq(Derivative(w1(t), t), -1.25662558596311w1(t))

= 2 阶方程 = dz/dt 方程: Eq(Derivative(z2(t), t), -0.176270821444344w1(t)**2 + 0.361648167900887w1(t)z1(t) - 0.0270616997342404z1(t)2 + 0.636625585963109z2(t)) dw/dt 方程: Eq(Derivative(w2(t), t), -0.34793830026576w1(t)2 + 0.713851832099113w1(t)z1(t) - 1.25662558596311w2(t) - 0.0534166785556559z1(t)**2)

= 3 阶方程 = dz/dt 方程: Eq(Derivative(z3(t), t), 0.0822183310251369w1(t)**3 - 0.309553239051093w1(t)2z1(t) - 0.352541642888688w1(t)w2(t) + 0.261022442368344w1(t)*z1(t)2 + 0.361648167900887w1(t)z2(t) + 0.361648167900887w2(t)z1(t) + 0.0585438170471337z1(t)**3 - 0.0541233994684808z1(t)z2(t) + 0.636625585963109z3(t)) dw/dt 方程: Eq(Derivative(w3(t), t), 0.1622895162862w1(t)**3 - 0.611022442368345w1(t)2z1(t) - 0.695876600531519w1(t)w2(t) + 0.515228239051092w1(t)*z1(t)2 + 0.713851832099113w1(t)z2(t) + 0.713851832099113w2(t)z1(t) - 1.25662558596311w3(t) + 0.115558752308197z1(t)*3 - 0.106833357111312z1(t)*z2(t))

= 4 阶方程 = dz/dt 方程: Eq(Derivative(z4(t), t), 0.01523321322074w1(t)**4 - 0.0849983517112508w1(t)3z1(t) + 0.246654993075411w1(t)2w2(t) + 0.147219050217673w1(t)2*z1(t)2 - 0.309553239051093w1(t)**2z2(t) - 0.619106478102185w1(t)w2(t)z1(t) - 0.352541642888688w1(t)w3(t) - 0.0562855430568364w1(t)z1(t)**3 + 0.522044884736689w1(t)z1(t)z2(t) + 0.361648167900887w1(t)z3(t) - 0.176270821444344w2(t)**2 + 0.261022442368344w2(t)z1(t)**2 + 0.361648167900887w2(t)z2(t) + 0.361648167900887w3(t)z1(t) - 0.0400490772251903z1(t)4 + 0.175631451141401*z1(t)2z2(t) - 0.0541233994684808z1(t)z3(t) - 0.0270616997342404z2(t)2 + 0.636625585963109z4(t)) dw/dt 方程: Eq(Derivative(w4(t), t), 0.0300686084751904w1(t)4 - 0.167776956943164w1(t)**3z1(t) + 0.486868548858599w1(t)**2w2(t) + 0.290593449782327w1(t)**2z1(t)2 - 0.611022442368345*w1(t)2z2(t) - 1.22204488473669w1(t)w2(t)z1(t) - 0.695876600531519w1(t)w3(t) - 0.111101179538749w1(t)z1(t)3 + 1.03045647810218w1(t)z1(t)z2(t) + 0.713851832099113w1(t)z3(t) - 0.34793830026576w2(t)2 + 0.515228239051092w2(t)z1(t)2 + 0.713851832099113w2(t)z2(t) + 0.713851832099113w3(t)z1(t) - 1.25662558596311w4(t) - 0.079052265955115z1(t)4 + 0.346676256924589z1(t)**2z2(t) - 0.106833357111312z1(t)z3(t) - 0.053416678555656*z2(t)**2)

= 5 阶方程 = dz/dt 方程: Eq(Derivative(z5(t), t), -0.0038619196646498w1(t)**5 + 0.028757566314031w1(t)4z1(t) + 0.06093285288296w1(t)3w2(t) - 0.0765875581962067w1(t)3*z1(t)2 - 0.0849983517112508w1(t)**3z2(t) - 0.254995055133752w1(t)**2w2(t)z1(t) + 0.246654993075411w1(t)2w3(t) + 0.0782583525562788w1(t)2z1(t)**3 + 0.294438100435346w1(t)2z1(t)z2(t) - 0.309553239051093*w1(t)2z3(t) + 0.246654993075411w1(t)w2(t)**2 + 0.294438100435346w1(t)w2(t)z1(t)2 - 0.619106478102185w1(t)w2(t)z2(t) - 0.619106478102185w1(t)w3(t)z1(t) - 0.352541642888688w1(t)w4(t) - 0.00527193636971858w1(t)z1(t)4 - 0.168856629170509w1(t)z1(t)2z2(t) + 0.522044884736689w1(t)z1(t)z3(t) + 0.261022442368344w1(t)z2(t)2 + 0.361648167900887w1(t)z4(t) - 0.309553239051093w2(t)**2z1(t) - 0.352541642888688w2(t)w3(t) - 0.0562855430568364w2(t)z1(t)3 + 0.522044884736689w2(t)z1(t)z2(t) + 0.361648167900887w2(t)z3(t) + 0.261022442368344w3(t)*z1(t)2 + 0.361648167900887w3(t)z2(t) + 0.361648167900887w4(t)z1(t) - 0.0263287080709955z1(t)**5 - 0.160196308900761z1(t)3z2(t) + 0.175631451141401z1(t)2z3(t) + 0.175631451141401z1(t)z2(t)**2 - 0.0541233994684807z1(t)z4(t) - 0.0541233994684807z2(t)z3(t) + 0.636625585963109z5(t)) dw/dt 方程: Eq(Derivative(w5(t), t), -0.00762298463733787w1(t)**5 + 0.0567641238697187w1(t)4z1(t) + 0.120274433900761w1(t)3w2(t) - 0.151175019222945w1(t)3*z1(t)2 - 0.167776956943164w1(t)**3z2(t) - 0.503330870829491w1(t)**2w2(t)z1(t) + 0.486868548858599w1(t)2w3(t) + 0.154472974862873w1(t)2z1(t)**3 + 0.581186899564655w1(t)2z1(t)z2(t) - 0.611022442368344*w1(t)2z3(t) + 0.486868548858599w1(t)w2(t)**2 + 0.581186899564655w1(t)w2(t)z1(t)2 - 1.22204488473669w1(t)w2(t)z2(t) - 1.22204488473669w1(t)w3(t)z1(t) - 0.69587660053152w1(t)w4(t) - 0.0104061952202809w1(t)z1(t)4 - 0.333303538616248w1(t)z1(t)2z2(t) + 1.03045647810218w1(t)z1(t)z3(t) + 0.515228239051092w1(t)z2(t)2 + 0.713851832099113w1(t)z4(t) - 0.611022442368344w2(t)**2z1(t) - 0.69587660053152w2(t)w3(t) - 0.111101179538749w2(t)z1(t)3 + 1.03045647810218w2(t)z1(t)z2(t) + 0.713851832099113w2(t)z3(t) + 0.515228239051092w3(t)*z1(t)2 + 0.713851832099113w3(t)z2(t) + 0.713851832099113w4(t)z1(t) - 1.25662558596311w5(t) - 0.0519698374317043z1(t)5 - 0.31620906382046*z1(t)3z2(t) + 0.346676256924589z1(t)2z3(t) + 0.346676256924589z1(t)*z2(t)2 - 0.106833357111312z1(t)z4(t) - 0.106833357111312z2(t)z3(t))

求解各阶解，形如：

= 各阶解 = w1 = 0 w2 = -0.0211143402100825a(t0)**2exp(1.27325117192622t)exp(-1.27325117192622t0) w3 = -0.00179505802924157a(t0)3exp(1.27325117192622t)exp(-1.27325117192622t0) + 0.0331683056091246*a(t0)3exp(1.90987675788933t)exp(-1.90987675788933t0) w4 = 0.00157395461517941a(t0)**4exp(1.27325117192622t)exp(-1.27325117192622t0) + 0.00422975802518084a(t0)4exp(1.90987675788933t)exp(-1.90987675788933t0) - 0.0223672954802892*a(t0)4exp(2.54650234385244t)exp(-2.54650234385244t0) w5 = -0.000525489787814235a(t0)**5exp(1.27325117192622t)exp(-1.27325117192622t0) - 0.00361886497519084a(t0)5exp(1.90987675788933t)exp(-1.90987675788933t0) - 0.00380315870113136*a(t0)5exp(2.54650234385244t)exp(-2.54650234385244t0) - 0.00401881299356255a(t0)**5exp(3.18312792981555t)exp(-3.18312792981555t0) z1 = a(t0)exp(0.636625585963109t)exp(-0.636625585963109t0) z2 = 0.0425080303571217a(t0)2exp(0.636625585963109t)exp(-0.636625585963109t0) - 0.0425080303571217*a(t0)2exp(1.27325117192622t)exp(-1.27325117192622t0) z3 = -0.0381756377357578a(t0)**3exp(0.636625585963109t)exp(-0.636625585963109t0) - 0.00361386528968398a(t0)3exp(1.27325117192622t)exp(-1.27325117192622t0) + 0.0417895030254418*a(t0)3exp(1.90987675788933t)exp(-1.90987675788933t0) z4 = 0.0140666785450764a(t0)**4exp(0.636625585963109t)exp(-0.636625585963109t0) + 0.00316873318782794a(t0)4exp(1.27325117192622t)exp(-1.27325117192622t0) + 0.00532916838964353*a(t0)4exp(1.90987675788933t)exp(-1.90987675788933t0) - 0.0225645801225478a(t0)**4exp(2.54650234385244t)exp(-2.54650234385244t0) z5 = 0.0130038185344948a(t0)5exp(0.636625585963109t)exp(-0.636625585963109t0) - 0.00105793198511116*a(t0)5exp(1.27325117192622t)exp(-1.27325117192622t0) - 0.00455949033428469a(t0)**5exp(1.90987675788933t)exp(-1.90987675788933t0) - 0.00383670342737988a(t0)5exp(2.54650234385244t)exp(-2.54650234385244t0) - 0.00354969278771912*a(t0)5exp(3.18312792981555t)exp(-3.18312792981555t0)

3.重整化方程

整合 $z, w$ 的所有解集，对于 $z (t)$ ，当 $t \to t_{0}$ ，计算 $\frac{d a (t_{0})}{d t_{0}}$ ，得到重整化方程，形如：

\begin{matrix} \frac{d a (t_{0})}{d t_{0}} = - 0.0228458146881948 a^{5} - 0.0342929006076182 a^{4} + 0.0509078545937214 a^{3} \\ - 0.0270616997342403 a^{2} + 0.636625585963109 a \end{matrix}

至此，原方程被转换为一维变尺度高阶动力学方程，对原方程有很好的近似，可以更加精确地刻画吸引域边界

故障检测

如图所示，红色曲线为故障曲线，当其突破吸引域边界（黑色点线）后一去不返，此时恢复系统参数无法使系统回归正常状态（stable point），利用重整化近似得到的边界（蓝色）能够及时检测到截断时间 $c c t$ ，而主流的线性近似方法（绿色虚线，IEEE Trans提出）则无法精确检测到。